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Abstract 

The Debye charging method is generalized to study the linear response properties of the asym¬ 
metric primitive model for electrolytes. Analytic results are obtained for the effective charge 
distributions of constituent ions inside the electrolyte, from which all static linear response prop¬ 
erties of system follow. It is found that, as the ion density increases, both the screening length and 
the dielectric constant receive substantial renormalization due to ionic correlations. Furthermore, 
the valence of larger ion is substantially renormalized upwards by ionic correlations, whilst that 
of smaller ions remains approximately the same. For sufficiently high density, the system exhibit 
charge oscillations. The threshold ion density for charge oscillation is much lower than the cor¬ 
responding value for symmetric electrolytes. Our results agree well with large scale Monte Carlo 
simulations. 

PACS numbers: 
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I. INTRODUCTION 


It was pointed out by Kirkwood [I| long ago that in sufficiently high densities, the ion- 
ion correlation functions of a symmetric electrolyte decays in an oscillatory fashion, a phe¬ 
nomenon called “charge oscillation” or “charge ordering”. Both analytic and numerical 
methods have been applied to study this phenomenon. At the threshold of oscillation, the 
Debye length comparable with the ion diameter, and hence it was often argued that the 
mechanism of charge oscillation is the competition between hardcore repulsion and Coulom- 
bic attraction between opposite ions. One-component plasma (OOP) also exhibits charge 
oscillation at high density [2|. The underlying mechanism is the strong electrostatic re¬ 
pulsion between (likely-charged) ions, much like in dense neutral liquids. Since OOP can 
be understood as the limit of extremely asymmetric electrolyte, where the valence of one 
component goes to zero, whilst symmetric electrolytes can be understood as the symmetric 
limit of asymmetric electrolytes, we would expect that charge oscillation also appears in 
asymmetric electrolytes, and the underlying mechanism is a combination of hard repulsion 
and the electrostatic repulsion between higher valence ions. 

Theoretically, the essence of charge oscillation can be captured by the renormalized elec¬ 
trostatic Green’s function G R (x — y), which is defined as the mean potential at x, due to a 
unit charge fixed at y and all other screening ions. The far held behaviors of all ion-ion pair 
correlations functions (which are more frequently used in liquid state physics) are identical 
to those of G R (x — y). The Green’s function however has the simplicity of satisfying a linear 
equation in the whole space: 

(-A + a*)G R {x~y) = ^S(x-y). (1.1) 

This equation was first derived by Kjellander and Mitchell Pffiu, in the setting of “dressed- 
ion theory”. The kernel a(x — y) can be expressed in terms of charge-charge correlation 
functions [6]. In this work, however, we shall not need this relation. As shown by Kjellander 
and Mitchell, the Green’s function decays in the form of a screened Coulomb potential in 
the far field: 

e -K R \x-y\ 

^e R \x-y\' 

The parameters k r and e R are determined by the pole structure of the Fourier transform 
a(k), and are generically different from the dielectric constant of the pure solvent and the 


( 1 . 2 ) 


G R (x — y) 
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bare inverse Debye length as given in PB theory. In the charge oscillation regime, rr and Cr 


become complex valued, and Eq. (1.2) should be understood as taking the real part. 


More remarkably, Kjellander and Mitchell have also shown that the mean potential d, 


due to a fixed constituent ion of specie /i satisfies an equation similar to Eq. (1.1), but with 
a renormalized charge distribution: 


(-A + a*) <fc(z - y) = -K^x - y), 


whose solution has the following far field asymptotics: 

q^e~ KRr 


(1.3) 


,(3-y) = 


4ne, 


(1.4) 


R 


where plays the role of renormalized charge of the ion. The physical significance of K^(f) 
is the effective charge distribution of ions of specie /i. There is an exact relation between the 


kernels a and all Kf s: 


0 


a = — 
e 




(1.5a) 


where q^n^ are, respectively, the bare charge and bulk density of ions of specie /j. The 


Fourier transforms of a, K., are related to Rr, €r and q^f via 


k r — o.{±iR R ), 


1 


2 i£rR-r 
r 


= Res 


Gn(k),iKR 


= K^IRr). 


(1.5b) 

(1.5c) 

(1.5d) 


Setting k = irr in the Fourier transform of Eq. (1.5a), we find 

0 


k r ~ 




R 

ifi • 


( 1 . 6 ) 


In the Poisson-Boltzmann theory, we have K M = q^S(r), i.e., the effective charge distribution 
of an ion is just the bare charge distribution. Hence = q^, Er = £, and Rr = Rq, where 

0 


r 2 0 = - 




(1.7) 


Eqs. (1.1 1.7) summarize the main results of “dressed-ion theory” due to Kjellander and 


Mitchell. According to these relations, all linear response properties of the electrolyte, in¬ 
cluding renormalizations of charges, Debye length, as well as dielectric constant are encoded 


3 












in the set of effective charge distributions K tl . For reviews of the dressed ion theory, see 
on H 2 ]. Note that our notations are different from those of Kjellander and Mitchell. In 
reference |E], this theory is reformulated in a form that can describe ion-specific interactions. 

The main purpose of this work is to compute approximately the effective charge dis¬ 
tributions K jl: in asymmetric primitive model. We shall develop an analytic formalism for 
effective charge distribution of a generic charged hard sphere particle immersed in an elec¬ 
trolyte (Sec. 0. By identifying this particle with a constituent ion, we obtain the effective 
charge distribution K fl for each specie of ions, and further the linear response kernel a, as 
well as various other parameters, e.g. qjf, kr, €r. As a special case, we shall apply the formal¬ 
ism to symmetric electrolytes, compute various renormalized parameters, and compare with 


previous theoretical results (Sec. III). We shall also apply the formalism to asymmetric elec¬ 


trolyte (Sec. IV). In both cases, we compare our analytic results with large scale numerical 
simulations and find good agreements. 


II. ANALYTIC FORMALISM 


A. The Method of Debye Charging 


We shall study the primitive model of electrolytes, where the solvent is modeled implicitly 
as a homogeneous media with a dielectric constant e, whilst ions are modeled as hard spheres 
with the same diameter d, and with a point charge at the center. There is no non-electrostatic 
interactions other than the volume exclusion. Furthermore, we assume that there is one 
specie of positive ion with charges q + = me and one specie of negative ion with charge 
g_ = — ne. Hence m,n are the valences of the positive and negative ions, whereas e is the 
fundamental unit of electric charge. Condition of charge neutrality then requires 

m p + — n p_ = 0. (2.1) 


The Hamiltonian of a homogeneous unperturbed electrolyte is 


Ho 


^ G o(Xij) + fHC i\Xij\/d) , 

i<j 


( 2 . 2 ) 


where g* and x % are the charge and position of i-th constituent ion, = x 3 — x 3 is the 
relative coordinate between q t , q 3 , and Go(r) = l/Aner is the electrostatic Green’s function 
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FIG. 1: Insertion of a hard sphere (blue disk) with diameter D creates a spherical excluded region 
for the centers-of-mass of all mobile ions, with a radius R c = (D + d)/ 2. The surface of the excluded 
region is called the surface of contact , schematically illustrated as the dashed circle. 

in the bulk solvent, whilst the function Uhc(£) describes the hardcore interaction: 

f oo, 0 < £ < 1; 

fHc(0 = < (2.3) 

l 0, £ > 1. 

The canonical partition function and the Helmholtz free energy of the homogeneous elec¬ 
trolyte are given by 

N 

Z 0 = Tre~ pH °= / e ~^ H °, (2.4a) 

J i= 1 

F 0 = —k B T log Zq = -k B T log Tr e~^ H °. (2.4b) 

We shall perturb the homogeneous electrolyte by two means simultaneously: 1) Intro¬ 
ducing an external potential £> ex , which is sufficiently weak so that it can be treated using 
the linear response theory; 2) Inserting at fa hard sphere particle with diameter D and 
with a point charge Q located at the center. The center-of-mass coordinates of all mobile 
ions are consequently constrained outside a sphere with radius R c = (D + d)/ 2. This sphere 
shall be called the contact surface, and R c its radius. For an illustration, see Fig. [TJ The 
total Hamiltonian of the perturbed system is then give by: 

H = H 0 + J2 x (^) + Q 0 ex (F) + Hep. (2.5) 

i 

The second and third terms in r.li.s. are, respectively, the interaction between the external 
potential 0 ex and the electrolyte, and that between 0 ex and the inserted particle. The last 
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term is the interaction between the inserted particle and all mobile ions: 

Hep = ^2 ^QqiG 0 (xi - r) + n H c(|^i - r\/R c ). 

i 

The change of free energy due to insertion of the particle is 

A F(r,Q,R c )= - i B T logTre _e ("" +E * ,rf " w+w " w +«.sp) 


( 2 . 6 ) 


+ ksT log Tr e 


-/3(-ff 0 +JV <?A ex Gd) 


= - k B T log 


e -0(£i ?^ ex (^)+Q^- ex (f)+^Ep) \ 

_/o 

= -/3(E i<li<t> e *(xi)) 


(2.7) 


where (• )o means average over the Gibbs distribution e /3 ^° of the unperturbed electrolyte. 
We shall call A F(f, Q, R c ) the free energy of insertion. 


In conjunction with the model system Eq. (2.5), we shall also consider two other related 
systems: 


1. Electrolyte perturbed by </> ex , but with no particle inserted: 


h = h 0 + J2 


qrf [Xi 


( 2 . 8 ) 


All mobile ions in the electrolyte respond to the external potential, and create a total 
mean potential 4>(r), which, according to the dressed-ion theory, satisfies a linear 


integro-differential equation [c.f. Eq. (1.1) 


- A0(f) + a * 4>(r) = 


(2.9) 


where p ex (r) = — eA</> ex (r) is the external charge distribution that generates 0 ex (r) in 
the first place. 


2. Electrolyte perturbed by the test particle, but with 0 ex (r) switched off: 


H = H 0 + H e p . 


( 2 . 10 ) 


The corresponding free energy of insertion can be obtained from Eq. (2.7) by setting 
0 ex (r) = 0: 


lim^A F(f,Q,R c ) -k B T log(e i3HEP ) 0 - 


( 2 . 11 ) 


Note that it is independent of the location of insertion. 
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The difference between Eqs. (2.7) and (2.11) is the potential of mean force (PMF) of the 
test particle, i.e., the effective interaction between the particle Q and the external potential 
(ft ex (r)\ 


U(f,Q,R c ) = A F{f,Q,R c )~ lim A F(r,Q,R c ) 

</> ex —>0 


= -k B T log 


e -j8(£i 5i0 ex (^)+Q</> ex (O+^ E p) \ 

/o 


( 2 . 12 ) 

(2.13) 


If the potential 0 ex (r) vanishes in the bulk (which is necessarily correct, if the external charge 
distribution p ex (r) is localized), U(r,Q]<f>) is also the work needed to bring the particle from 
the bulk to the present position x. This is in fact the definition of PMF used in most 
literatures. 

For a weak potential, U(x, Q, R c ) can be expressed as a linear functional of (ft ex . Actually, 
it is more useful to express U(x,Q,R c ) in terms of the total mean potential (ft prior to the 
insertion (which satisfies Eq. (0): 

U(x,Q,R c )= iKix-iJ.qmy,. (2.14) 

■'y 

Note that (ft and (ft ex are proportional to each other in the linear regime. The kernel K(x — 
y, Q) represents the effective charge density of the inserted particle. For the simple case of 
an infinitesimal point charge, the PMF is U(x,dq, 0) = dq(ft(x), which corresponds to an 
effective charge distribution K(x — y) = dq 5{x — y). 

To simplify the analysis, we shall choose the mean potential to be monochromatic <ft(r) = 


Ak-r 


(22! Additionally, we shall also choose to insert the particle at the origin, so Eq. (|2.14|) 


becomes 


U(0,Q,R c ) = K(k,Q)(fto. 


(2.15) 


Hence, all we have to do is to calculate the PMF [7(0, Q, R c ) of the test particle in a 
monochromatic background potential. 


Let us now take the derivative of Eq. (2.7) with respect to Q: 

Tr [JA qiG 0 (xi - r) + (ft ex (r )] e~ pH 


^AF { f,Q,R c ) = 


Tre-Pn 

QiGoAi - + (ft ex (r) 

= AgQ), 


(2.16) 
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where the average (•) is defined with respect to the Gibbs measure e l3H , with H given by 


Eq. (2.5). Hence 0( f,Q ) as defined is the mean potential acting on Q at the center of the 


test particle, due to all mobile charges {qi} as well as the external charges p ex (f). Now if 


we switch off the external potential in Eq. (2.16) (and exchange the limit and derivative): 


d 

— lim A F(f,Q,R c )= lim 0(f,Q). 
dQ </> ex ^o 7 v ' 


(2.17) 


Subtracting this from Eq. (2.16) and using (2.12), we find: 


r\ r\ r\ 

U(r,Q,R c ) = AF(f,Q)-— lim AF(r,Q) 

= ip(r,Q)~ lim 0(f, Q) = 50(r, Q), 
<^> ex —^0 


(2.18) 


where 5if(r,Q) is the difference between and its bulk value. Following Debye [18]. 


we integrate Eq. (2.18) over Q, we find 


n Q 


U (r, Q, R c ) = U (r, 0, R c ) + / dqdifif.q). 

Jo 


(2.19a) 


Combining this result with Eq. (2.15), we see that K(k,Q) can be calculated if we know 
-0(0, Q) (the mean potential acting on the test particle at the origin) as well as t/(0, 0, R c ), 


i.e., the PMF of a neutral particle. We shall discuss 0(0, Q) in Sec. IIB and Sec. IIC and 


then discuss U(r,0,R c ) in Sec. IID 


B. Mean potential acting on the test particle 


To compute 0(0, Q), we shall first compute the total mean potential $(r, Q) at r, due to 
both the test ion fixed at the origin and all other mobile ions. 0(0, Q) can then be obtained 
from <3>(r, Q) by subtracting off the Coulomb potential due to Q itself and taking the local 
limit: 


0(0, Q) = lim 


r—>0 


mQ) 


Q 

47rer 


( 2 . 20 ) 


Inside the contact surface r < R Cl no other ions can enter, and hence the potential <l>(r', Q ) 
satisfies the Poisson equation: 

- eV 2< h(r, Q)=Q 5(f), r < R c . (2.21a) 


Outside the contact surface r > R c , we shall use Eq. (2.14) to approximate the PMF of all 


constituent ions[23], so that the ion number density of specie p is 

W0P) = fq, e' 


= ^ p -PK^{r,Q) 


(2.21b) 















Consequently, <L(r, Q) satisfies the following nonlinear (and nonlocal) partial integro- 
differential equation: 


- eV 2 $(f, Q) = Y, % + /f(f), r > R c 


(2.21c) 


Note that linearization of Eq. (2.21c) (together with Eq. (1.5a)) leads to Eq. (2.9). 


Eq. (2.21c) is an improvement over the nonlinear PBE, and reduces to the latter if one 
approximate K t , by Additionally, <f>(r, Q) satisfies the standard electrostatic bound¬ 

ary conditions: 


lim $(r, Q) = 0, 

r—» oo 

^ d 

$(r, Q ), —<h(r, Q ) continuous at r = R c . 
or 


(2.22a) 

(2.22b) 


We shall calculate the PMF up to the second order in Q. In view of Eq. (2.19), we only 


need to solve Eqs. (2.21c) to the Erst order in Q and in </> 0 . We decompose <L(r, Q ) into four 
parts: 

<D(f,Q) = (t>{r) + 0 h ' c '(f) + Q [G(r) + 0 c (r)] + 0(Q 2 ) + O(0 2 ), (2.23) 

where <fi(r) is the mean potential in the absence of the test ion, and satisfies the linear integro- 
differential equation (2.9). 0 hx '(r) arises due to the insertion of a neutral hard sphere. G(r) 
is independent of </>, whilst (j) c is linear in q i. Both G{r) and (j) c are independent of Q. All 
the ignored terms are at least quadratic either in Q or in (j). 


Let us set Q = 0 in Eq. (2.23), and obtain 

<P(r, 0) = 0(r) + 0 h ' c '(r) + O(0 2 ). 


(2.24) 


This satisfies Eqs. (2.21a) and ( 2.21c[ ) with Q = 0, and corresponds to inserting a neutral 
hard sphere at the origin. Expanding these equations to first order in 0 and 0 hx ', and 


subtracting off Eq. (2.9), we find: 

- A0 hx ' 


r = 


—a * 0 hx -(r), 
a * (ft(r), 


r > R c ; 
r < R r . 


(2.25a) 


The equation satisfied by G(r) can be obtained by substituting Eq. (2.23) into Eqs. (2.21a) 


and (2.21c) and extract the part that is linear in Q and independent of (j): 

-AG(r) + a * G(r) = 0, r > R c ; 


-AG(r) = -5(f), 
e 


(2.25b) 


r < R r . 
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One must be careful not identifying G(r) with the renormalized Green’s function G R (r) 


which satisfies Eq. (1.1). As one can see, G(r) explicitly take into account the effects of 


hardcore repulsion, whilst G R {f) does not. 

Finally the equation satisfied by 0 C can be obtained by extracting the bilinear term 


(proportional to Q4>) of Eqs. (2.21a) and (2.21c): 


-A0 c (r) + a * 0 c (f) = e l f3 2 ^ 

* [0(f) + 0 hx -(f)j A' t *G(r), r>R c ; 

—A0 c (r) = 0, r < R c . 

The boundary conditions for 0 hx '(r), G(r), and 0 c (f) are identical to those for $(F, Q) 


(2.25c) 


C. Local approximation 


Eqs. (2.25) are difficult to solve, because of the non-local nature of convolutions appearing 


in them. To simplify the problem, we shall make the following local approximation for the 
kernel a : 

a(k) = a(x) = n r5(x). (2.26a) 


Correspondingly, the Green’s function in Eqs. (1.1) are approximated by 

1 


G R (k) 


k 2 + k r ’ 


G R (f) 


47rer 




(2.26b) 


In another word, this amounts to approximate the renormalized Green’s function by a 
screened Coulomb potential with a renormalized Debye length. Such an approximation 
is motivated by two considerations: 1) to preserve the large scale feature of renormalized 
theory outside the hard core, and 2) to make the analytical calculation feasible. One can in 
principle make a more refined approximation, at a cost that the analyses can no longer be 


carried out explicitly. To preserve the exact relation between a and A0, Eq. (1.5a), we must 
make a similar approximation to the effective charge distributions A),: 


K^k) « Qr(<a), A0(f) « q R (q^)S(f). 


(2.26c) 


This amounts to simply replacing the bare charges by the renormalized charges, and ignoring 
the diffusive nature of the effective charge distributions. Such an approximation turns out 
be rather successful, as we shall demonstrate below. 
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With the local approximation, Eq. (2.25a) now reduces to 


A 0 1 


h.c. / 


-«r 0 hh r>R c ; 


r = 


40(r), 


(2.27) 


r < ify. 


Recall that 0(r) = 0o e* fc ' r has the form of plane wave, and can be expanded in terms of 
spherical harmonics Y) m using the well-known formula: 

OO l 


gik-r _ 


47T EE i l ji(kr)Y] m (k)Y lm (f) 

1=0 m=—l 

= jo(kr) + anistropic, 


(2.28) 


where ji(kr) are the spherical Bessel functions of the first kind, and k,f are the unit vectors 
parallel to k, r respectively. Likewise, 0 hx '(r) can also be expanded in terms of spherical 


harmonics. To satisfy Eq. (2.27) and be compatible with the boundary conditions at the 


origin and at the infinity, the expansion must have the following forms: 


/h.c., 


r > R c \ 


r = 


J2i, m a i k i{ K R r ) Y im(k)Y lm (r), 

Yji,m [°1 rl + M r )] Yim{k)Y lm (r), r < R c , 


(2.29) 


where ki(hiRr) are the modified spherical Bessel functions of the second type, which vanish 
as v — y oo , and ai,ci,bi(r ) must be found. Because 0 h,c -(r) is continuous at the origin, 
all functions bi(r) with nonzero l must vanish at r = 0. m We shall need two pieces 
of information about 0 hx ’(E) in this work: 1) 0 hx, (O) and 2) 0 hx, (r) averaged over the 
contact surface r = R c . For both quantities, the anisotropic channels with l fy 0 make no 


contribution. For the isotropic channel, l = 0 and Eq. (2.27) reduces to 


f d? .2 d | ,h.c./ \ , ,2 /h.c. 


1 d^ + rd~r'^ ' C ' (r) + KR(t) X ' (r) = °’ r > Rc] 




(2.30) 


r < R r , 


The solution is (after imposing the boundary conditions Eq. (2.22)) 


L h.c., 


r = 


{ 0o fi(kR c , RrR c ) k 0 (R R r), 

k 2 

00 f 2 (kR c , RrRc) + -j§jo(kr) 


r > R c ; 


r < R c , 


(2.31) 
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where the functions fi(x,y),f 2 (x,y) are given by 


h(x,y) = 

f2(x,y) = 


y 3 e y (sin x — x cos x) 
x 3 (l + y) 

y 2 (y sin x + xcosx) 
£ 3 (1 + y) 


Using the same approximation Eqs. (|2.26|), Eq. (J2.25b[) can be easily solved: 

r > R c ; 
r < R r , 


G(r) = < 


47t(1 + K R R c )r ’ 


+ 


kr 


k 47rer 47re(l + krR c ) ’ 


Finally Eq. (2.25c) reduces to 




= e 


1 p 2 '52paQa (<&? [<l>(r) + 0 hx '(r)] G(r), r > R c ; 

r < R c . 


, d 2 2d. 

~ 1 *3 + M (r) - °. 


(2.32a) 

(2.32b) 


(2.33) 


(2.34) 


Eq. (2.34) will be solved in Sec. IV for asymmetric electrolytes. For symmetry electrolytes, 


</> c (r) vanishes identically due to symmetry reason. 


D. PMF of a neutral particle: “contact value theorem” 

We shall now outline a method for the PMF of a neutral hard sphere inside an electrolyte, 
Eq. (2.15). As illustrated in Fig. [lj all ions are excluded from the region r < R c . Hence the 
partition function of the total system is 


Z = 


e~ 0H = 


'Rc 


jWn^e 


'Y[d 3 fid(ri-R c ) 

where fj is the unit vector parallel to Tj. The total free energy is: 

F = —ksT log Z = Fq + AF(0, 0, R c ). 


2 ~. P ~BH 


(2.35) 


(2.36) 


where F 0 is the free energy of the unperturbed electrolyte, whilst AF(0,0,i? c ) is the free 


energy of insertion of the neutral hard sphere (c.f. Eq. (2.7)). Now let us take the differential 
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of F with respect to R c : 


OF 

OR, 


dR c = dRc^~Yl fd 2 '^ II fd 3 Rj 0( r j — R ) (e pH ) r .- 

i 

= k B T dR c <f)d 2 f/ <5(r) — r) 

i 

= k B T dR c j) d 2 f ^ n /t (r), 


(2.37) 


where n^if) is the average ion number density of species /q and the integral j> d 2 r is over the 


contact surface. Eq. (2.37) is a variation of the contact value theorem |T7], which gives an 
exact relation between the particle number density for hard sphere systems and the pressure 


acting on a hard wall. It is important to note in Eq. (2.37), we have used the fact that 


the Hamiltonian is independent of R c . This would not be correct if there is image charge 
effects. Luckily enough, in the primitive model, the dielectric constants of the ions and of 
the solvent are identical, and hence image charge interactions do not appear. 


Now, according to Eq. (2.21b), the ion number density n^{f) (with a neutral particle 
fixed at the origin) is given by: 




(2.38) 


where $(r, 0) is given in Eq. (2.24). Substituting this back into Eq. (2.37), linearizing in 


terms of <L(r, 0), using the local approximation Eq. (2.26c), and further integrating over 


the contact surface, we can express the r.li.s of Eq. (2.37) as a linear functional of <f>(r, 0) 


Finally keeping the part that is linear in <L(r, 0), and integrating over the radius of contact 
surface R Cl we obtain the PMF of a neutral hard sphere: 

-Rc 


U (0,0, R c ) = —An V'' [ 

Jo 


dR c R 2 c {$(f, 0)) cont , 


(2.39) 


where (• ) con t means average over the contact surface. We re-emphasize that this result is 
applicable only if the hard sphere has the same dielectric constant as the solvent. 


III. SYMMETRIC ELECTROLYTES 

Let us first apply the general method to the simple case of symmetric electrolytes, where 
positive/negatives ions have charges ±g and hard sphere diameter d. The renormalized 
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charges of positive and negative ions remain opposite to each other: — Qr(— q) = qn{q) = Or- 


Consequently the r.h.s. of (2.34) vanishes identically. This means 0 C vanishes identically, to 


the first order in 0. Likewise, the r.h.s. of Eq. (2.39) vanishes identically. Hence {7(0, 0, R c ) = 

0. 


The total potential $(r, Q) can be obtained by substituting Eqs. (2.31), (2.33) back into 


Eq. (2.23). Using Eq. (2.20), we further calculate 0(0 ,Q), the mean potential acting on Q: 


"0(0, Q) = </>o 




1 + + f2{kR c , KrRc) 


Q K-R 


4ire(l+K R R c 


(3.1) 


We now use Eq. (2.18) to compute 50(0, Q), and use Eq. (2.19) and the fact that K(k, 0) 


vanishes to obtain the effective charge distribution K(k,Q ) for the test particle: 

m Q) =q {i+1 - 4 }+°m 


(3.2) 


which is an entire function of k. Recall that we have calculated K(k,Q) up to the order 
of Q 2 , hence the ignored terms are at least of order Q 3 . The renormalized charge can be 


obtained using Eq. (1.5d): 


Qr(Q) — K (■ iKR , Q) — 


Qe 


KrRc 


0(Q 3 ). 


(3.3) 


1 + KrR c 

The fact that there is no term of order of Q 2 is actually enforced by the charge inversion 
symmetry: Qr(-Q) = -Qr(Q). 


A. Renormalized Debye length and renormalized dielectric constant 


We may apply Eq. (3.2) to the constituent ions with R c = d,Q = q, and further apply 


Eq. (1.5a) to calculate the linear response kernel a: 


a(k) = til 1 + t§ + f 2 (kd, K R d) 


(3.4) 


where f 2 (x,y) was already defined in Eqs. (2.32) 


The renormalized Debye length can be obtained via Eqs. (1.5) and (3.3): 

2 


« o 


„K R d 


Or 


(3.5a) 


1 + Krid q 

where qR is the renormalized charge of the positive constituent ion. This is a self-consistent 
equation for the renormalized inverse Debye length kr. We can also calculate the renormalize 
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dielectric constant: 


— = 2 - -K R d - e~ KRd 
e 2 


(1 + n R d) - - sinh n R d 


(3.5b) 


B. Charge Oscillation 


Careful analysis of Eq. ( |3.5a[ ) indicates a critical value Kq defined by 

1.3465, 


n'„d = \/2 (2 + V3) e i 'V ss 


(3.6) 


such that for k 0 > Kq, there is no real root for Eq. (3.5a). What happens is that a pair of 
real roots collide with each other at Kq an d bifurcate into the complex plane. Consequently 
the renormalized charge q R and renormalized dielectric constant e R also become complex 
valued. Whilst our local approximation is not really applicable if n 0 is sufficiently close 


to Kq: d seems very natural to argue that one should take the real parts of Eq.. (1.2) and 


(1.4) if the relevant parameters become complex. Therefore in the regime Ko > Kq, mean 
potential decays in an oscillatory fashion and the system exhibits charge oscillation. The 
corresponding Kj t at the threshold is 


K R d = 1 + y/3 


2.732. 


(3.7a) 


C. Comparison with simulations 

We simulated 1 : — 1 electrolytes with three different ion sizes: d = 5A, 7.5A, 10A re¬ 
spectively, and determine all renormalized parameters. The simulation method is described 
elsewhere. [16] As shown in Fig. [2j our MC results seem to agree with Eqs. (3.5) remark¬ 
ably well, both below and above threshold of charge oscillation. A general argument due 
to Kjcllander and Mitchell shows that e R should vanish at the threshold. This contradicts 


our result Eq. (3.5b), which gives ~ 0.64. We note that near the threshold, the local 


approximation Eq. (2.26b) breaks down, and therefore Eq. (3.5b) can not be trusted. In 


any case, simulations near the threshold are very difficult and we have no reliable results to 
report so far. 


15 


















K R d 



Im(A R d) 


(a) 



K 0 d 


e R 



K 0 d 


FIG. 2: Renormalized parameters for dense symmetric electrolyte. Left: Renormalized 
v.s. bare inverse Debye length. Middle: the imaginary part of krcI in the charge oscillation regime. 


Symbols: MC simulation results. Solid curves: our analytic results Eq. (3.3). The dashed curves 


are, respectively: black, generalized Debye-Huckel by Lee and Fisher [J>|; red and blue: results using 
other approaches, LMPB and MSA, both from reference [10j. The other popular theory HNC does 
not yield a close form result. Right: Renormalized dielectric constant, for which we have not found 
any previous analytic result. Simulations details are discussed in a separate publication [16] . 

D. Effective charge distribution and mean potential 

The real space version of the effective charge distribution can also be calculated, by 

,Kr(krR c ~ KR'r + 1 ) 


Fourier transforming Eq. (3.2): 


K (r, Q) = Q 5(r) + Q- 


-e(Rc-r), 


(3.8) 


4ir(l + k r R c )r 

where 9(R c — r) is the Heaviside step function. The first term (Dirac delta function) is clearly 
due to the bare charge of the test particle. The second term is positive and monotonically 
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decreasing, and vanishes identically outside the contact surface (r > R c ). It is the diffusive 
part due to charge correlations. Note that even though the second term is singular (diverges 
as r -1 ) at the origin, it does not generates any singularity in the mean potential. The mean 


potential, which is related to K(k,Q ) via Eq. (1.3), can also be explicitly calculated: 

Q R e- KRr 


<!>(?, Q) = 


Aner 

Q 


Qk r 


r > R c ] 


0 < r < R c , 


(3.9) 


47rer 47re(l + krR c ) ’ 

which has the exact form of screened Coulomb potential outside the contact surface. Since 


Qr is given by Eq. (3.3) and always has the same sign as the bare charge, we see that there 


is no charge inversion in symmetry electrolytes at the level of our approximation. 

The results shown in this section were also obtained by Kjcllander [2UJ some time ago. 
We rederive these results to illustrate the method of Debye charging. 


IV. ASYMMETRIC ELECTROLYTES 


Let us now consider asymmetric electrolytes. Let the renormalized charges of constituent 
ions be, respectively, q R = m R e, and q R = Ur, where m R = q+/e, n R = q R /e are the 
renormalized valences. Using the neutrality condition (2.1) in Eqs. E3 and ( |1.6[ ) we can 
obtain 


Kg = e 1 f3e 2 p + m(m + n), 
Kr = e~ 1 f3e 2 p + m(m R +n R ). 


Dividing, we find a useful relation between kr and the renormalized valences: 

2 


kr 

k 0 


m R + n R 
m + n 


(4.1a) 

(4.1b) 

(4.2) 


A. PMF of neutral hard sphere continued 


Let us calculate the PMF of a neutral hard sphere U(0,0,R c ) using Eq. (2.39). For 
this purpose, we need the mean potential $(r, 0) in the presence of a neutral hard sphere, 


Eq. (2.24), with 0 hx -(r) given by Eq. (2.29). Averaging over the contact surface is trivial, 
because we have already thrown out the anisotropic parts. The result is 


(<L(f, 0)) cont = 0o [jo(kRc) + fi{kR c , K R R c )k 0 (K R R c )]. 


(4.3) 
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Substituting this into Eq. (2.39) and further back into Eq. (2.15), we find the effective charge 


distribution K(k, 0) for a neutral hard sphere: 


K(k,0) = 

where the function \&o(x,i/) is defined as 


47ri? c ^5^ PaQa') 


ty()(kRc, KrR c ). 


(x,y) = / [j Q (xt) + fi(xt,yt) k 0 {yt)]t 2 dt 
J o 

= —j— — x Ci (x + xr/ _1 ) (x cos(x?/ _1 ) — y sin (an/ -1 )) 
x y i 

+ xCi ( xy _1 ) (xcos (xy -1 ) — y sin (xy -1 )) 

— x Si (x + xy -1 ) (x sin ( xy _1 ) + y cos (xy *)) 

+ xSi(xr/ _1 ) (x sin (xy x ) + y cos ( X V x )) 

— y ((x 2 + 2 y) cos(x) + x(y — 2) sin(x)) + 2 y 


where Ci(x) and Si(x) are the cosine integral and sine integral functions: 

Ci(s) = - 


* 00 cost f z sin t 

-dt, Si(z) = / —:— dt. 


t 


'0 


t 


(4.4) 


(4.5) 


(4.6) 


Both Ci(^) and Si(z), and hence ^o(kR c , krR c ) as well, are entire functions. 

Using Eqs. (2. 1|) and (4.1), we can rewrite Eq. (4.4) into the following dimensionless form: 


U(M) = (4.7) 

e b mnym + n) 

where b = e 2 /AneT is the Bjerrum length. If both the inserted particle and the constituent 
ions are point-like, R c —>■ 0, and To(0,0) = 1/3, and hence 


47T 


K(k,0) -> — 


(4.8) 


B. The correlation potential </> c (r) 


To calculate </> c (0), we only need to solve the isotopic channel of Eq. (2.34). The isotropic 


component of r.li.s. of Eq. (2.34) can be easily shown as 


x 


<hQ 0 , td \ r ) 2 , n c 

- 9{r - R c )(m R - n R )K R {K R b) — -— 

e 1 + krR, 

ko(i^Rr) [jo(kr) + fi(kR c , K R R c )k 0 (K R r )]. 


(4.9) 
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where fi(x,y) is defined in Eqs. (2.32), and jo(u) = sin u/u, ko(u) = e u /u. 0 c (r) can 
now be found using standard Liouvillc method 


0 C (O) = 


2 c f> 0 Q 


(m R - n R ) (K R b) ^ 2 {kR c , k r R c ), 


^2 (x,y) = 


=2 y 


Vi(x,y) Ei(3y) 


4(1 + y) 2 

+ i (-) (El(2 y + ix) - E x (2 y - ix )) 


x 


(4.10) 


(4.11) 


where is Ei(z) the exponential integral function EU, defined as 


Ei ( 2 ) = / t e dt. 


(4.12) 


Ei(^) has a logarithmic singularity at the origin z — 0, and a branch cut on the negative 
real axis. Consequently, the function ^ 2 (x,y) as a function of complex variable x has two 
branch cuts on the imaginary axis: one from 2 iy to ioo, and the other from —2 iy to —zoo. 
These singularities have no influence on the leading order asymptotics of mean potential by 
a charged hard sphere, or on the interaction between two charged spheres, as the leading 
order asymptotics of the latter quantities are controlled by the pole k = in R , where the 
function 'h 2 {kR c , krR c ) is analytic. 

Substituting this and Eqs. (2.31), ( |2.33[ ) back into Eqs. ( |2.23 ) and (2.20), we find the 
potential acting on the the charge Q (subtracting its bulk value): 


6rp(0,Q) = 0o 
2 Q 




1 + -gf + f 2 {kR c , krR c ) 


H- (rn R - n R ) ( K R b ) ^ 2 (kR c , k r R c ) 

e 


(4.13) 


Now using Eqs. (2.19) and (4.7) to carry out the Debye charging process, we finally obtain 
the effective charge distribution (in Fourier space): 


1 T> ( r K o ( k R R c) 3 ( m R n ~ n R m ) lTf /?t-) „ 

~K(k, Q) = - 2^ — n - 7 - ; — ^o[kR c ,K R R c 

e k r [RrO) mnym + n) 


(4.14) 


+ j [l + + fz{kd, n R d )] 

RJ \ 2 

—J {m R - n R ) ( K R b ) ^ 2 (kR c , r r R c ), 


where To, T 2 are defined in Eqs. (4.5) and (4.11). 
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C. Renormalized Debye length of asymmetric electrolytes 


We can now set Q = me, —ne in Eq. (4.14) to obtain the effective charge distributions 


for each specie of constituent ions in the bulk and use Eq. (1.5a) to find the linear response 


kernel a(k). We shall skip the calculation details and present the results directly. The kernel 
a(k) is 


a(k) = nl 1 + + f 2 (kR c , k r R c ) 

+ (m - n)(m R - n R ){K R b)^ 2 (kR c , k r R c ) 


(4.15) 


On the other hand, setting k = in R in Eq. (4.14), and using Eq. (1.5d), we End the renor¬ 


malized charge of a hard sphere with bare charge Q: 

Qr _ (Q\ (Q\ 2 

— — ao + cq I — 1 + a 2 I — 1 + 0(Q ). 

where the coefficients ao,ai,a 2 are 

k 2 0 1 (i m R n - n R m) 


(4.16a) 


o-o — 


x 


a i = 


Kft (R R b) mn(m + n ) 


—Ei (—k r R c — 1) H—Ei(—1) + e KRRc (n R R c — 2) + 2 
e e 

P^rRc 


1 + r r R c 
02 = ( m R - n R ) ( K R b)e 2KRR ' 
Ei (k r R„ 


(4.16b) 


x 


a (rrRc ~ l)e 2Kj?i?c Ei(3 k r R c 


[A(k r R c + l) 2 4 (k r R c + l) 3 

Note that for symmetric electrolytes, the even order coefficients ao, a 2 are contained to vanish 
by symmetry, since m R = n R , and m = n. The lowest order renormalization is therefore of 


order Q i , see Eq. (3.3) 


Setting Q = me,—ne, the l.h.s. of Eq. (4.16a) reduces to the renormalized valences 


m R , — n R of the constituent ions. Solving for m R , n R we hnd: 


rn R = (4.17a) 

F\{R R d)m [F 0 (hi R d)hil/+ mns{F 2 {n R d)sn{m + n) — 1)] 
mns [F 2 (n R d)s (m 2 + n 2 ) — 1] — [ F 2 (n R d)s{m — n) 2 — 1] F 0 (n R d)i% 1/n R : 
tir = (4.17b) 

Fi(K R d)n[F 0 (K R d)KQ/K R + mns(F 2 (K R d)sm(m + n) — 1)] 
mns [F 2 (K R d)s (m 2 + n 2 ) — 1] — [ F 2 {n R d)s{m — n) 2 — 1] F 0 (K R d)i% 1/tt R 


20 





















(e) (f) 

FIG. 3: Renormalized parameters for asymmetric electrolyte. Comparison of theoretical 


predictions Eqs. (4.17) with large scale MC simulations. Panels (a), (c), (e) : 2 : —1 ,d = 7.5 A. 
The straight-line in panel (a) is the prediction of classical PB. Panels (b), (d), (f): 3 : — 1, d = 2AA. 
Note that the threshold densities for charge oscillations for these systems are much lower than 
that for symmetric electrolytes. Substantial disagreements between theory and simulation for cr 
arise near the threshold of charge oscillation. This is probably due to the local approximation 


Eq. (2.26a). Simulations details are discussed in a separate publication US- 
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where s = krI) and functions F 0 (y), Fi(y), F 2 (y) are defined as 


F 0 {y) = [- Ei(—1 — y) — - Ei(l) + e y (2 — y) — 2 

y A le e 

e y 

My) = 


i + y 


(4.17c) 


My) = 


3 2 y 


4(1 + y) 3 L 


(1 - y) e 2y Ei(3y) + (1 + y)E 1 (y) 


Plugging these back into Eq. (4.2), we find a self-consistent equation for k r \ 

M 2 

K 0 


(4.17d) 


M K R d ) 


Fq(k,rcI) ((r() + mns [ 2 F 2 (nRd)mns — 1 


mns [ F 2 (KRd)s (m 2 + n 2 ) — l] — [F 2 (nRd)s(m — n ) 2 — l] F 0 (K R d) 

This gives the renormalized inverse Debye length kr as an implicit function of the bare one 


Kq. Finally the renormalized dielectric constant can also be obtained using Eq. (1.5) (with 
shorthand y = KRd below): 

e R_ l= i 

(4.17e) 


24 (y + l) 3 \k r J 
12 (y + l) 2 [ 2 + 2 y — 3e y + y sinh y + (y + 1) cosh y ] 

+ s e~ y (‘m — n)(m R — n R ) (y + 1) [3e 3y Ei(— y) + 3e 2y + 1 ] 


3e 3y Ei(—3y) [ 2e y ( y 2 + 3) sinh y — 6 ye v cosh y + y + 1 ] 


In Fig. [3] we show the comparison between the analytic results Eqs. (4.17) and large scale 
MC simulations of two systems: 1) 2 : —1 ,d = 7.5 A, and 2) 3 : —1 ,d = 24 A. It can be 
seen that the agreement between theory and simulations is generally good. Note that as the 
density increases, the valence of larger ion is substantially renormalized upwards by ionic 
correlations, whilst that of smaller ions remains approximately the same. Furthermore, 
both systems exhibit charge oscillations in the high density regime. The threshold value 
Kq d is approximately unity, not too much different from that of symmetric electrolytes. 


However, because the bare Debye length is related to ion valences nonlinearly via Eq. (1.7), 


the threshold ion density for charge oscillation is much lower in asymmetric electrolytes than 
in asymmetric electrolytes. 
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Let us now check the special case of symmetric electrolytes, where m = n, uir = Ur. All 


results reduce to those in Sec. Ill In particular, Eqs. (4.17a)-(4.17d) reduce to Eq. (3.5a), 


and Eq. (4.17e) reduces to Eq. (3.5b), as it should be. 

Let us also check the point-like ion limit, where d —>■ 0. The renormalized charge 


Eq. (4.16a) reduces to the following limiting form: 


Or _ Q log 3 

e e 4 


(n R b)(m R - n R ) I — ) + 0(Q 3 ). 


(4.18) 


Hence charge renormalization becomes significant if the renormalized Debye length becomes 
comparable with the Bjerrum length, i.e., when (/%&) ~ 1. In the same limit, Eqs. (4.17c) 


reduce to F 0 —> 0, Fi —> 1, F 2 —> jlog3, and hence Eqs. (4.17a) and (4.17b) reduce to 


Eq. (4.17d) reduces to 


m R 

4 — n(K,Rb)(m + n) log(3) 

(4.19) 

m 

4 — («r&) (m 2 + n 2 ) log(3) ’ 

riR 

4 — m(nRb)(m + n) log(3) 

(4.20) 

n 

4 — ( kr 5 ) (m 2 + n 2 ) log(3) ’ 

4 

4 — 2{KRb)mn\og{?>) 

(4.21) 

4 

4 — (krF) (m 2 + n 2 ) log(3) ’ 


and finally Eq. (4.17e) reduces to 


e_R _ _ (fifcfe) (4 - 3 log 3 ) (m - n) 2 

e 12(2 — mn(K R b ) log 3 ) 


(4.22) 


These results also predict charge oscillation if (kr&) is comparable with unity, which clearly 
indicate that charge oscillation in asymmetric electrolytes can be solely driven by electro¬ 
static correlations (mainly between ions of higher valences). 

Strictly speaking, a two-component plasma model with point-like ions is not well-defined 
because of the instability towards annihilation of opposite charges. For small but finite ion 
sizes, this instability is manifested as formation of bound pairs of ions, or even larger clusters. 
This instability does not show up at the level of our approximation, just as in the classical 
Debye-Huckel theory. Linearization, which is adopted in both theories, is responsible for the 
suppression of this instability at short scales. As a logical consequence, whenever bound ion 
clusters can not be ignored, linearization breaks down, and the short scale properties derived 
in our theory (similar to those of Debye-Huckel theory) can not be trusted. This happens, 
for example, for dense electrolytes with small ions where the electrostatic energy between 
neighboring ions becomes much larger than the thermal energy Indeed, we have 
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also simulated 3 : —1 electrolyte with d = 7.5 A. The largest interaction energy between 
two ions (in close contact) is then approximately SkgT, which makes linearization a bad 
approximation in the dense regime. As expected, we found that substantial disagreement 
between theory and simulation. 
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Our results demonstrate that renormalization of ion valences and Debye length, dielectric 
constant is more apparent in asymmetric electrolytes than in symmetric. We also find that, 
generically, the valence of larger ions is renormalized substantially upwards, whereas that of 
smaller ions remains stays approximately constant. Finally, the threshold density for charge 
oscillation in asymmetric electrolytes is much lower than that for symmetric electrolytes. 
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